/******************************************************************************
 * 
 * Announce: CSharpKit, Basic algorithms, components and definitions.
 *           Copyright (C) ShenYongchen.
 *           All rights reserved.
 *   Author: 申永辰.郑州 (shenyczz@163.com)
 *  WebSite: http://gitee.com/shenyczz/csharp_csk
 *
 * THIS CODE IS LICENSED UNDER THE MIT LICENSE (MIT).
 * THIS CODE IS PROVIDED *AS IS* WITHOUT WARRANTY OF 
 * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING ANY
 * IMPLIED WARRANTIES OF FITNESS FOR A PARTICULAR
 * PURPOSE, MERCHANTABILITY, OR NON-INFRINGEMENT.
 * 
******************************************************************************/

using System;

namespace CSharpKit.Numerics.Utils
{
    /// <summary>
    /// Barnes 带通滤波 (高斯加权客观分析)
    /// </summary>
    public class KBarnes
    {
        public KBarnes()
        {
            // 参数选取C1=10000，C2=150000，G=0.3主要保留了400-1000公里的波长的中尺度系统。
        }

        public int Rows { get; set; }
        public int Cols { get; set; }

        double x0 { get; set; }
        double y0 { get; set; }
        double xInterval { get; set; }
        double yInterval { get; set; }

        // 滤波参数
        /// <summary>
        /// 滤波常数，通常取0.2-0.4
        /// </summary>
        /// <value></value>
        double G { get; set; }
        /// <summary>
        /// 
        /// </summary>
        /// <value></value>
        double C { get; set; }

        double[] _X { get; set; }
        double[] _Y { get; set; }


        double[] Datas { get; set; }
        double[] m_ppD;			// 误差

        //  0-第一次滤波场 1-第二次滤波场
        double[] m_ppF0;			// 第一次滤波结果
        double[] m_ppF1;			// 第二次滤波结果

        /// <summary>
        /// 设置格点参数
        /// </summary>
        /// <param name="datas">X方向格点数</param>
        /// <param name="nx">X方向格点数</param>
        /// <param name="ny">Y方向格点数</param>
        /// <param name="x0">网格区域左下角经度</param>
        /// <param name="y0">网格区域左下角纬度</param>
        /// <param name="xInterval">经度网格间距</param>
        /// <param name="yInterval">纬度网格间距</param>
        public void SetGrid(double[] datas, int nx, int ny, double x0, double y0, double xInterval, double yInterval)
        {
            this.Datas = datas;
            this.Cols = nx;           //列
            this.Rows = ny;           //行

            this.x0 = x0;          //网格区域左下角经度
            this.y0 = y0;          //网格区域左下角纬度

            this.xInterval = xInterval;    //经度网格间距
            this.yInterval = yInterval;    //纬度网格间距

            this.m_ppD = new double[nx * ny];
            this.m_ppF0 = new double[nx * ny];
            this.m_ppF1 = new double[nx * ny];
        }

        void SetParam_Filter(float g, float c)
        {
            G = g; C = c;
        }


        /// <summary>
        /// c1,g1,c2,g2
        /// </summary>
        public void Transact()
        {
            // 第一次滤波,计算初值场
            for (int i = 0; i < Rows; i++)
            {
                var y = y0 + i * yInterval;
                for (int j = 0; j < Cols; j++)
                {
                    var x = x0 + j * xInterval;
                    // 
                    var pos = i * Cols + j;
                    m_ppF0[pos] = (float)Calculate(x, y, 1.0, C, Datas);	//f0
                    m_ppD[pos] = Datas[pos] - m_ppF0[pos];			//误差
                }
            }

            for (int i = 0; i < Rows; i++)
            {
                var y = y0 + i * yInterval;
                for (int j = 0; j < Cols; j++)
                {
                    var x = x0 + j * xInterval;
                    // 
                    var pos = i * Cols + j;
                    m_ppF1[pos] = (float)Calculate(x, y, G, C, m_ppD); //插值误差
                    m_ppF1[pos] = m_ppF0[pos] + m_ppF1[pos];			//f1
                }
            }

        }


        /// <summary>
        /// 
        /// </summary>
        /// <param name="x0"></param>
        /// <param name="y0"></param>
        /// <param name="g">滤波参数 g</param>
        /// <param name="c">滤波参数 c</param>
        /// <param name="ppData"></param>
        /// <returns></returns>
        double Calculate(double x0, double y0, double g, double c, double[] ppData)
        {
            double sum_w = 0.0;
            double sum_wfk = 0.0;

            for (int i = 0; i < Rows; i++)
            {
                var y = y0 + i * yInterval;
                for (int j = 0; j < Cols; j++)
                {
                    var x = x0 + j * xInterval;

                    var pos = i * Cols + j;
                    var fk = ppData[pos];		                // 原始场数据
                    var r = _distance(x0, y0, x, y) * 111.0;   //计算距离,一个经纬度为111.0km
                    var w = _gauss(r, g, c);	                //计算权重系数

                    sum_w += w; // 权重系数和
                    sum_wfk += fk * w;		//原始场数据权重值的和
                }
            }

            return (sum_wfk / sum_w);
        }


        double _distance(double x0, double y0, double x1, double y1)
        {
            return System.Math.Sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0));
        }

        /// <summary>
        /// 高斯权重
        /// </summary>
        /// <param name="r">距离</param>
        /// <param name="g"></param>
        /// <param name="c"></param>
        /// <returns></returns>
        double _gauss(double r, double g, double c)
        {
            return System.Math.Exp(-r * r / (4.0 * g * c));
        }


        // {{END}}
    }

}
